str(cat$countries)

pays<-cat$countries[cat$countries$income_level!="Aggregates",c("iso3c", "country","capital_city","longitude","latitude", "region","income_level")]
kable(head(pays))

indic<-cat$indicators
kable(head(indic,3))

vars <- wbsearch(pattern = "CO2",fields="indicator")
kable(vars)

meta<-cat$indicators[cat$indicators$indicator_id %in% c("EN.ATM.CO2E.KT","EN.ATM.CO2E.PC"),]
kable(meta)

tabCO2 <- wb_data(indicator = c("EN.ATM.CO2E.KT"),
                 return_wide = TRUE,
                 mrv=1,
                 country ="countries_only")



tabCO2<-tabCO2[order(tabCO2$EN.ATM.CO2E.KT,decreasing = T),]
kable(head(tabCO2,10))

tabCO2 <- wb_data(indicator = c("EN.ATM.CO2E.PC"),
                 return_wide = TRUE,
                 mrv=1,
                 country ="countries_only")



tabCO2<-tabCO2[order(tabCO2$EN.ATM.CO2E.PC,decreasing = T),]
kable(head(tabCO2,10))

world_data <- wb_data(indicator = c("AG.SRF.TOTL.K2","AG.LND.ARBL.HA","SP.POP.TOTL","SP.URB.TOTL","NY.GDP.MKTP.CD", "EN.ATM.CO2E.KT"),
                 return_wide = TRUE,
                 start_date = 1990,
                 end_date = 2019,
                 country ="countries_only")

world_data<-world_data[,-c(1,3)]

# recodage (attention : ordre alphabetique)
names(world_data)<-c("iso3c","date","ARB","SRF","CO2","GDP","POP","URB")

pays<-cat$countries[cat$countries$income_level!="Aggregates",c("iso3c", "country","capital_city","longitude","latitude", "region","income_level")]
tab<-merge(pays,world_data, by="iso3c",all.x=F,all.y=T)
kable(head(tab))

saveRDS(tab,"data/wb_don_1990_2019.Rdata")

# Extract meta
meta<-cat$indicators[cat$indicators$indicator_id %in% c("AG.SRF.TOTL.K2","AG.LND.ARBL.HA","SP.POP.TOTL","SP.URB.TOTL","NY.GDP.MKTP.CD", "EN.ATM.CO2E.KT"),]

# Select information and add personal code
meta<-data.frame(meta[,c(1,2,4,5)])
meta$shortcode<-c("ARB","SRF","CO2","GDP","POP","URB")
meta<-meta[,c(5,1,2,3,4)]

# Display
kable(meta)

# Save meta
saveRDS(meta,"data/wb_don_1990_2019_meta.Rdata")
write.table(meta,"data/wb_don_1990_2019_meta.csv",sep=";",dec=",", row.names = F,fileEncoding = "UTF-8")

# Load data
don<-readRDS("data/wb_don_1990_2019.Rdata")



tabmis<-don %>% select(iso3c,country,date,ARB,SRF,CO2,GDP,POP,URB)  %>%
  gather("ARB", "SRF","CO2","GDP","POP","URB",key = VAR, value = DON) %>%
  mutate(Missing = is.na(DON))

kable(head(tabmis))             

res <- tabmis %>% group_by(VAR, date) %>% 
                  summarise(nbmis=sum(Missing), nb=n(), pct = 100*nbmis/nb) 

p<-ggplot(res) + 
  aes(x=date,y=pct,color=VAR) +
  geom_line() +
  scale_y_continuous("% de pays à valeurs manquantes", breaks = c(0,10,20,30,40,50,60,70,80,90,100))+
  scale_x_continuous("Année", breaks = c(1990, 1995, 2000, 2005, 2010, 2015, 2020))
p

res <- tabmis %>% group_by(iso3c, country) %>% 
                  summarise(nbmis=sum(Missing), nb=n(), pct = 100*nbmis/nb) %>%
                  arrange(-pct)

kable(head(res,20))

# load a data.frame with columns space, time, indic 
estim<-function(df   = don)               # dataframe
 
{


# Select variable
sel<-dcast(df,space~time)
M<-as.matrix(sel[,-1])
rownames(M)<-sel$space

####################################################
########   Estimation Procedure ####################
####################################################

# check for lines or column  filled with missing values and eliminate
check1<-apply(!is.na(M),FUN="sum",1)
check2<-apply(!is.na(M),FUN="sum",2)
M<-M[check1>0,check2>0]



dim(M)


# define dimensions
nr<-nrow(M)
nr
nc<-ncol(M)
nc


# create time matrix with NA
Mt<-matrix(rep(1:nc,nr),nrow=nr,nc=nc,byrow=T)
rownames(Mt)<-rownames(M)
colnames(Mt)<-colnames(M)
Mt[is.na(M)]<-NA
Mt[1:4,1:4]



# identify previous and next available valuer
Mt_min<-Mt
Mt_max<-Mt
for ( i in 1:nr) { 
  for ( j in 1:nc) { 
    if (is.na(Mt[i,j])) {
      Mt_min[i,j]<-max(Mt[i,1:j],na.rm=T)
      Mt_max[i,j]<-min(Mt[i,j:nc],na.rm=T)    
    }  
  }
}
Mt_min[1:4,1:4]
Mt_max[1:4,1:4]


# choose estimation method
M_met<-matrix("OK",nrow=nr,ncol=nc)
rownames(M_met)<-rownames(M)
colnames(M_met)<-colnames(M)
for ( i in 1:nr) { 
  for ( j in 1:nc) { 
    if (is.na(Mt[i,j])) {
      M_met[i,j]<-"IN"
      if (is.infinite(Mt_max[i,j])) {M_met[i,j]<-"EX"} 
      if (is.infinite(Mt_min[i,j])) {M_met[i,j]<-"RE"}  
    }   
  }
}
M_met[1:4,1:4]

###### Estimation of missing values  ####

M_est<-M
str(M_est)

#### step 1: Interpolation ######

for ( i in 1:nr) { 
  for ( j in 1:nc) { 
    if ((M_met[i,j]=="IN")) {
      t0<-Mt_min[i,j]
      t1<-Mt_max[i,j] 
      tacm<-(M[i,t1]/M[i,t0])**(1/(t1-t0))
      M_est[i,j]<-M[i,t0]*(tacm**(j-t0))
    }   
  }
}
M_est[1:4,1:4]


#### step 2: Extrapolation ######

for ( i in 1:nr) { 
  for ( j in 2:nc) { 
    if ((M_met[i,j]=="EX")) {
      N<-M_est
      t0<-j-1
      t1<-j
      N<-M_est[is.na(M_est[,t0])==F,]
      N<-N[is.na(N[,t1])==F,]
      tacm<-(sum(N[,t1])/sum(N[,t0]))**(1/(t1-t0)) 
      M_est[i,j]<-M_est[i,t0]*(tacm**(t1-t0))
    }   
  }
}
M_est[1:4,1:4]


#### step 3: Retropolation ######

for ( i in 1:nr) { 
  for ( j in (nc-1):1) { 
    if ((M_met[i,j]=="RE")) {
      N<-M_est
      t0<-j
      t1<-j+1
      N<-M_est[is.na(M_est[,t0])==F,]
      N<-N[is.na(N[,t1])==F,]
      tacm<-(sum(N[,t1])/sum(N[,t0]))**(1/(t1-t0)) 
      M_est[i,j]<-M_est[i,t1]/(tacm**(t1-t0))
    }   
  }
}
M_est[1:4,1:4]



#######################################################
############# EXPORT RESULTS  #########################
#######################################################

x<-reshape2::melt(M_est)
names(x)<-c("space","time","estim")
y<-reshape2::melt(M_met)
names(y)<-c("space","time","method")
z<-merge(x,y, by=c("space","time"))

return(z)
}

# Load data
don<-readRDS("data/wb_don_1990_2019.Rdata")

# Estim ARB
df<-data.table(space=don$iso3c,time=don$date,index=don$ARB)
est<-estim(df)[,1:3]

names(est)<-c("iso3c","date","ARB_est")
don<-merge(don,est, by=c("iso3c","date"), all.x=T,al.y=F)

# Estim SRF
df<-data.table(space=don$iso3c,time=don$date,index=don$SRF)
est<-estim(df)[,1:3]

names(est)<-c("iso3c","date","SRF_est")
don<-merge(don,est, by=c("iso3c","date"), all.x=T,al.y=F)

# Estim CO2
df<-data.table(space=don$iso3c,time=don$date,index=don$CO2)
est<-estim(df)[,1:3]

names(est)<-c("iso3c","date","CO2_est")
don<-merge(don,est, by=c("iso3c","date"), all.x=T,al.y=F)

# Estim GDP
df<-data.table(space=don$iso3c,time=don$date,index=don$GDP)
est<-estim(df)[,1:3]

names(est)<-c("iso3c","date","GDP_est")
don<-merge(don,est, by=c("iso3c","date"), all.x=T,al.y=F)

# Estim POP
df<-data.table(space=don$iso3c,time=don$date,index=don$POP)
est<-estim(df)[,1:3]

names(est)<-c("iso3c","date","POP_est")
don<-merge(don,est, by=c("iso3c","date"), all.x=T,al.y=F)

# Estim URB
df<-data.table(space=don$iso3c,time=don$date,index=don$URB)
est<-estim(df)[,1:3]

names(est)<-c("iso3c","date","URB_est")
don<-merge(don,est, by=c("iso3c","date"), all.x=T,al.y=F)

# Save
saveRDS(don,"data/wb_don_1990_2019.Rdata")

don<-readRDS("data/wb_don_1990_2019.Rdata")
#don<-as_tibble(don)
don2<-don %>%  group_by(date) %>% mutate(ARB_rnk = rank(desc(ARB_est)),
                                         SRF_rnk = rank(desc(SRF_est)),
                                         CO2_rnk = rank(desc(CO2_est)),
                                         GDP_rnk = rank(desc(GDP_est)),
                                         POP_rnk = rank(desc(POP_est)),
                                         URB_rnk = rank(desc(URB_est)),
                                         TOT_rnk = (ARB_rnk+SRF_rnk+CO2_rnk+GDP_rnk+POP_rnk+URB_rnk)/6) 

tab<-don2 %>% filter(date==1996) %>% 
             select(iso3c, country, TOT_rnk, ARB_rnk,SRF_rnk,CO2_rnk,GDP_rnk,POP_rnk,URB_rnk,TOT_rnk) %>%
             arrange(TOT_rnk)
kable(head(tab,10), digits=1)

tab<-don2 %>% filter(date==2016) %>% 
             select(iso3c, country, TOT_rnk, ARB_rnk,SRF_rnk,CO2_rnk,GDP_rnk,POP_rnk,URB_rnk,TOT_rnk) %>%
             arrange(TOT_rnk)
kable(head(tab,10), digits=1)

don2<-don2 %>%  group_by(date) %>% mutate(ARB_wgt = 100*ARB_est/sum(ARB_est,na.rm=T),
                                         SRF_wgt = 100*SRF_est/sum(SRF_est,na.rm=T),
                                         CO2_wgt = 100*CO2_est/sum(CO2_est,na.rm=T),
                                         GDP_wgt = 100*GDP_est/sum(GDP_est,na.rm=T),
                                         POP_wgt = 100*POP_est/sum(POP_est,na.rm=T),
                                         URB_wgt = 100*URB_est/sum(URB_est,na.rm=T),
                                         TOT_wgt = (ARB_wgt+SRF_wgt+CO2_wgt+GDP_wgt+POP_wgt+URB_wgt)/6)

tab<-don2 %>% filter(date==1996) %>% 
             select(iso3c, country, TOT_wgt, ARB_wgt,SRF_wgt,CO2_wgt,GDP_wgt,POP_wgt,URB_wgt,TOT_wgt) %>%
             arrange(desc(TOT_wgt))
kable(head(tab,10), digits=1)

tab<-don2 %>% filter(date==2016) %>% 
             select(iso3c, country, TOT_wgt, ARB_wgt,SRF_wgt,CO2_wgt,GDP_wgt,POP_wgt,URB_wgt,TOT_wgt) %>%
             arrange(desc(TOT_wgt))
kable(head(tab,10), digits=1)

saveRDS(don2,"data/wb_don_1990_2019.Rdata")
write.table(don2,"data/wb_don_1990_2019.csv",sep=";",dec=",", row.names = F,fileEncoding = "UTF-8")

map<-st_as_sf(sovereignty110)
par(mar=c(0,0,0,0), mfrow=c(1,1))
plot(map$geometry,col="lightyellow")

dim(map)

sel<-map[map$adm0_a3 %in% c("FRA", "NCL","ATA","ATF","USA", "PRI","CHN","TWN","MAR", "SAH","CHN","TWN","ISR","PSX"),c("sovereignt","sov_a3","type","admin", "adm0_a3","name","note_adm0","iso_a3","wb_a3")]
kable(sel)

map<-st_as_sf(sovereignty110)
par(mar=c(0,0,0,0))
plot(map$geometry,col="lightyellow")

dim(map)

par(mfrow=c(1,2))
plot(map[map$sov_a3 %in% c("ISR","JOR","SYR","LBN","EGY"),]$geometry, col=c("gray80","orange","gray80","gray80","gray80"))
title("Limits of Israël",cex=0.5)
plot(map[map$sov_a3 %in% c("MAR","SAH","DZA","MRT"),]$geometry, col=c("gray70","orange","gray70","lightyellow"))
title("Limits of Morocco")

map<-st_as_sf(countries110)
small<-st_as_sf(tiny_countries110)
par(mar=c(0,0,0,0), mfrow=c(1,1))
plot(map$geometry,col="lightyellow")
plot(small$geometry,col="red", add=T)

map<-st_as_sf(countries50)

par(mar=c(0,0,0,0))
plot(map$geometry,col="lightyellow")

# Load map
map<-st_as_sf(countries110)
map<-map[c("adm0_a3","name")]
names(map)<-c("iso3c","name","geometry")

# Add polygons center
coo<-st_coordinates(st_centroid(map,of_largest_polygon = T))
map$Lon<-coo[,1]
map$Lat<-coo[,2]


# adjust some codes
map$iso3c[map$iso3c=="KOS"]<-"XKX"    # Kosovo
map$iso3c[map$iso3c=="PSX"]<-"PSE"    # Palestinian territories
map$iso3c[map$iso3c=="SDS"]<-"SSD"    # South Sudan

# Save
st_write(map,"data/wb_map_low.shp",delete_dsn=T)

saveRDS(map,"data/wb_map_low.Rdata")

# Load map
map<-st_as_sf(countries50)
map<-map[c("adm0_a3","name")]
names(map)<-c("iso3c","name","geometry")

# Add polygons center
coo<-st_coordinates(st_centroid(map,of_largest_polygon = T))
map$Lon<-coo[,1]
map$Lat<-coo[,2]


# adjust some codes
map$iso3c[map$iso3c=="KOS"]<-"XKX"    # Kosovo
map$iso3c[map$iso3c=="PSX"]<-"PSE"    # Palestinian territories
map$iso3c[map$iso3c=="SDS"]<-"SSD"    # South Sudan

# Save
st_write(map,"data/wb_map_high.shp",delete_dsn=T)

saveRDS(map,"data/wb_map_high.Rdata")

don<-readRDS("data/wb_don_1990_2019.Rdata")
tab<-don %>% filter(date > 1994, date < 2016) %>%
  select (iso3c, country, date, CO2_est, POP_est, CO2_rnk, CO2_wgt) %>% 
  filter(CO2_rnk <11) 

tab<-as.data.frame(tab)


p <- plot_ly(tab,
             x = ~date,
             y = ~CO2_wgt,
             color= ~country) %>%
             add_lines()%>%
         layout(title = "Top 10 des émissions de CO2 (1995-2015)",
         yaxis = list(title = "% des émissions mondiales de CO2", type="log"), 
         xaxis = list(title = "Année"))
p

don<-readRDS("data/wb_don_1990_2019.Rdata")
tab<-don %>% filter(date > 1994, date < 2016) %>%
  select (iso3c, country, date, GDP_est, POP_est, GDP_rnk, GDP_wgt) %>% 
  filter(GDP_rnk <21) 

tab<-as.data.frame(tab)


p <- plot_ly(tab,
             x = ~date,
             y = ~GDP_wgt,
             color= ~country) %>%
             add_lines()%>%
         layout(title = "Top 20 du PIB en ppa (1995-2015)",
         yaxis = list(title = "% du PIB mondial en ppa", type="log"), 
         xaxis = list(title = "Année"))
p

# Load map
map<-readRDS("data/wb_map_high.Rdata")

# Load don for 2015
don<-readRDS("data/wb_don_1990_2019.Rdata")
don<-don[don$date==2015,]

# Merge with map
mapdon<-merge(don,map,by="iso3c",all.x=F,all.y=T)

# Create map



g <- list(showframe = TRUE,
          framecolor= toRGB("gray50"),
          coastlinecolor = toRGB("black"),
          showland = TRUE,
          landcolor = toRGB("lightyellow"),
          showcountries = TRUE,
          countrycolor = toRGB("black"),
          countrywidth = 0.2,
 #        projection = list(type = 'azimuthal equal area'))
          projection = list(type = 'Mercator'))



p<- plot_geo(mapdon)%>%
  add_markers(x = ~Lon,
              y = ~Lat,
              sizes = c(0, 1000),
              size = ~CO2_wgt,
              color= "red",
              hoverinfo = "text",
              text = ~paste('Pays: ', country,
                            '<br /> % mondial : ',round(CO2_wgt,3),
                            '<br /> rang mondial : ',round(CO2_rnk,0))) %>%
  layout(geo = g,
         title = "Part des émissions mondiales de CO2 en 2015",
         width=800,
         height = 400)
p

# Load map
map<-readRDS("data/wb_map_high.Rdata")

# Load don for 2015
don<-readRDS("data/wb_don_1990_2019.Rdata")
don<-don[don$date==2015,]

# Merge with map
mapdon<-merge(don,map,by="iso3c",all.x=F,all.y=T)

# Create map



g <- list(showframe = TRUE,
          framecolor= toRGB("gray50"),
          coastlinecolor = toRGB("black"),
          showland = TRUE,
          landcolor = toRGB("lightyellow"),
          showcountries = TRUE,
          countrycolor = toRGB("black"),
          countrywidth = 0.2,
 #        projection = list(type = 'azimuthal equal area'))
          projection = list(type = 'Mercator'))



p<- plot_geo(mapdon)%>%
  add_markers(x = ~Lon,
              y = ~Lat,
              sizes = c(0, 1000),
              size = ~GDP_wgt,
              color= "red",
              hoverinfo = "text",
              text = ~paste('Pays: ', country,
                            '<br /> % mondial : ',round(GDP_wgt,3),
                            '<br /> rang mondial : ',round(GDP_rnk,0))) %>%
  layout(geo = g,
         title = "Part du PIB mondial (ppa) en 2015",
         width=800,
         height = 400)
p